Hydrologic Design 4501

Chapter 7: Unit Hydrograpgh and Routing

Ardeshir Ebtehaj
University of Minnesota
Table of Contents

1- Linear Reservoir and Convolution

A linear reservoir is a simple and meaningful representation of the watershed response to precipitation events. Continuity equation for a linear reservoir is written as follows:
Figure 1: Schematic of a linear reservoir. A reservoir is linear if we assumed that the storage function is linearly related to the outflow through the constant K.
where is the outlfow, is the inflow (rainfall excess), is change in the basin storage over time, and K represents a time scale that captures the residence time of the reservoir. Here, the storage is the sum the rainfall excess minus the outflow up to time t.
Now upon substitution, we obtain and have to solve a first-order ordinary differential equation (ODE):
Multiplying both sides by , one can get:
Now because of the product rule, we can rewrite the above equation as follows:
,
and thus, we have,
Therefore, we can obtain the outflow of a linear reservior as follows:

1-1 Concept of Convolution

The above solution states that the outflow at any time is a linear sum of the contributions by an initial flow (first term) that goes to zero over time and the integral term (second term) that represents the effects of the reservoir resident time Kon the outflow. The integral term is a convolution integral. Learning the concept of convolution is crucial for understanding the concept of the unit hydrograph.
Mathematical convolution for two generic functions and is defined as follows:
convolution1.gifConvolution2.gif
Figure 2: Visualization of two examples of the convolution integral. We can see that the convolution of a function with a kernel may result in translation (shift) and dilution of . We will show later on the the rainfall excess can be represented by and the kernal function or the unit hydrograph is a characteristic of the watershed. The result of the convolution of excess rainfall with the unit hydrograph will the be watershed outflow.

1-2 Instantaneous Unit Hydrograph (IUH)

Now looking back on the integral in the linear reservoir outflow equation, we can see that the inflow resembles and the convolution kernal resembles
and .
If we define the input rainfall excess as a unit impulse function as follows:
Then, we can define the impulse response function of a linear reservoir as follows:
The function is often called the instantaneous unit hydrograph (IUH), which is the response of a linear reservoir to a unit excess rainfall with an infinitesimally small duration.

1-3 Unit Hydrograph

Similar to the above formalism for a linear reservoir, a watershed has a response function to a unit pulse of excess precipitation --- albeit a more complex one. This response function is known as the unit hydrograph, which is a time-invariant characteristic of a watershed. More specifically, the unit hydrograph of a watershed is defined as tthe direct runoff hydrograph (DRH) resulting from a unit pulse of excess rainfall, which is defined as a unit depth of excess rainfall (i.e., 1 , 1 , etc) distributed uniformly over the entire watershed for an effective duration of .
Figure 3: The response of a linear reservoir to a unit impulse (left) resembles the response of a watershed, called unit hydrograph, to a unit pulse of rainfall excess (right), which is a unit amount of rainfall, uniformly distributed over the watershed for duration of .
Note that convolution is a linear operator
which means that we can superimpose the impacts of two pulses of excess rainfall over a watershed through linear combination of their individual outflow responses.
Figure 4: Superimposition of the response functions (i.e., re-scaled and shifted unit hydrographs) to excess rainfall hyetographs. For example, the above schematic shows the response to two pulses of excess rainfall, where .
Therefore, having the unit hydrograph, one can obtain the outflow of a watershed as convolution of excess rainfall with the unit hydrograph as follows:
Below is an example of hydrograph superposition used to obtain the outflow of a watershed in response to a 24-hr storm -- when a 6-hourly unit hydrograph is available.
Figure 5: Example of 6-hr Unit Hydrograph convolution (Credit: The COMET Program).

2- Computation of the Unit Hydrograph

Figure 6: Convolution of excess rainfall hyetograph (ERH) at intervals with the unit hydrograph of a basin, sampled at points. The convolution results in a Direct Runoff hydrograph (DRH) at points, where .
Since, hydrologic data is not continuous, we must use a discretized version of the convolution equation as follows:
For the example shown in Figure 6, we have
where , . Therefore, the above convolution can be expressed as a linear system of equation as follows:
Through convolution, one can obtain the above set of equations through flipping the rainfall excess hyetograph and slide it over the coordinates of the unit hydrograph and then compute the sum over multiplications.

2-1 Forward substitution and least-squares

We can solve the above set of equations through forward substitution. In other words, one can obtain and then use the second equation to obtain and so forth.
One can arrange all data of the above linear system of equations in a matrix form as follows:
and solve it for values to obtain the unit hydrograph from the observed direct runoff hydrograph of a precipitation excess event.
We can see that in the above linear system of equations we have more equations that unknown. Therefore, calculation of the unit hydrograph is often an overdetermined problem.
The above system of linear equation, which is resulted from a discrete convolution operation can be rewritten in matrix form as follows:
where is an column vector, is an matrix, and is an column vector.
For the provided example, the matrix form is as follows:
The matrix form can be easily implemented into a software library for numerical linear algebra. As previously explained, the equation is an over-determined linear system, since we have more equations than unknowns. Two of the main methods used to solve this system of equations for are:
1. Forward/backward substitution: Write out the first or last L equations. Solve the first or last equations and sequentially substitute and solve the remaining equations. In this formalism only L equations are used, which result in multiple solutions to the problem.
2. Least squares: Using linear algebra, we can derive a formula to calculate the least squares estimate of the unit hydrograph ordinates as follows:
where is the matrix transposition operator, which flips the elements of a matrix over its diagonal elements. Note that by least squares, we mean that the solution has minimum error variance amoung all other possible solutions and the solution is unique.
The above matrix formulation of the linear system can be solved efficiently using a technical programing language such as MATLAB. Note that, in the above formulation, we multiply both sides with to obtain a square matrix , which might be invertible. We need to note that is often a rectangular matrix with more rows than columns () and cannot be inverted.
The matrix is called pseudo inverse of . In MATLAB, computation of the unit hydrograph can be simply implelmented using either of the following commands: u = P\q or u = pinv(P)*q.
--------------------------------------------------------------------------------------------
Example Problem 7.1: Find the half-hour unit hydrograph using the excess rainfall hyetograph and direct runoff hydrograoh shown in the following table.
Solution: Given , and since thus .
% Least squares solution for Unit Hydrograph
clear
format short g
close all
t = 0:0.5:5;
% Direct runoff hydrograph (DRH)
q = [423 1923 5297 9131 10625 7834 3921 1846 1402 830 313]';
% Excess rainfall Hyetograph (ERH)
p = [1.06 1.93 1.81 0 0 0 0 0 0 0 0];
r = [1.06 0 0 0 0 0 0 0 0];
P = toeplitz(p,r);
% Unit Hydrograph
% u = inv(P'*P)*P'*q;
% or
u = pinv(P)*q;
%% Plots
yyaxis left
plot(t',q,'DisplayName','DRH','LineWidth',2)
hold on
plot(t(1:length(u)),u,'DisplayName','UH','LineWidth',2)
xlabel('\bf Time (hr)')
ylabel('\bf Discharge (cfs)')
yyaxis right
bar(t,p,'FaceAlpha',0.3,'DisplayName','Precipitation')
ylabel('\bf Precipitation excess(in)')
legend('show')
--------------------------------------------------------------------------------------------
Example Problem 7.2: Calculate the streamflow hydrograoh for a storm of 6 inch excess rainfall over three hours with hourly rate of p = [2,3,1] using the half hour unit hydrograph computed in the previous example and assume the base below is constant at 500 [cfs] throughout the flood. Check that the total depth of direct runoff is equal to the the total depth of excess precipitation for the watershed area of 7.03 miles square.
Solution:
%computation of the strearmflow hydrograph
clear
format short g
close all
dt = 0.5; % [hr]
t = 0:dt:5;
A = 7.03*5280^2; %[ft^2]
Q_base = 500; % [cfs]
u = [404, 1079,2343,2506,1460,453,381,274,173]; % [cfs/in]
p = [2,3,1]; % [in]
Q_direct = conv([0 0 u 0 0],p,'valid'); % [cfs]
Q = Q_direct+Q_base; % [cfs]
V_d = sum(Q_direct)*dt; % [cfs.hr]
V_d = V_d*3600/1; % [ft^3]
r_d = V_d/A*12; % [in]
Q_direct_1 = conv([0 0 u 0 0],[2 0 0],'valid');
Q_direct_2 = conv([0 0 u 0 0],[0 3 0],'valid');
Q_direct_3 = conv([0 0 u 0 0],[0 0 1],'valid');
yyaxis left
plot(t,repmat(Q_base,size(t)),'-ok','DisplayName','Base Flow','LineWidth',2)
hold on
plot(t,Q_base+Q_direct_1,'-or','DisplayName','Response to p = [2,0,0] [in]','LineWidth',2);
plot(t,Q_base+Q_direct_1+Q_direct_2,'-og','DisplayName','Response to p = [2,3,0] [in]','LineWidth',2);
plot(t,Q_base+Q_direct_1+Q_direct_2+Q_direct_3,'-ob','DisplayName','Response to p = [2,3,1] [in]','LineWidth',2);
xlabel('\bf Time (hr)')
ylabel('\bf Discharge (cfs)')
yyaxis right
bar(t,[p, zeros(1,8)],'FaceAlpha',0.3,'DisplayName','Precipitation')
ylabel('\bf Precipitation excess (in)')
legend('show')

2-2 S-Hydrograph

As we mentioned, the UH is a time invariant characteristic of a watershed in response to a unit pulse of excess rainfall with a specific duration . However, since the UH is obtained through a linear process, we can change a unit hydrograph from one duration to another. To that end, we need to derive the S-hydrograph based on the principle of superposition.
By definition, the S-hydrograph is the direct runoff hydrograph of a watershed resulting from a continuous excess rainfall at a constant rate of 1 [cm/hr] or [in/hr] for an infinite period of time.
In other words, the S-hydrograph is the response of the watershed to summation of an infinite number of unit pulse of excess rainfall with intensity and duration -- or a unit rainfall excess over an infinite period of time (Figure 7).
Figure 7: Conceptual representation of the S-hydrograph.
As a result, the unit hydrograph for duration, with intiger n values, can be obtained as follows:
Figure 8: Derivation of the unit hydrograph with duration , where , from an S-hydrograph that is generated from a unit hydrograph with duration .
--------------------------------------------------------------------------------------------
Example Problem 7.3: Use 0.5-hr unit hydrograph in the following table to calculate 1.5-hr unit hydrograph.
Solution:
A shifted version of the hydrogrsph is shown between columns 2 and 3. Column 3 represents the sum of the unit hydrographs multiplied by [hr]. Column 4 shows the lagged unit hydrograph for [hr]
An excel file solving an example to derive a 1.5-h unit hydrograph from 0.5-hr unit hydrograph is posted in the course website. Column 5 computes the unit hydrograph for [hr].
% computing 1.5-hr hydrograph from 0.5-hr hydrograph using S-hydrograph
format short g
clear
dt = 0.5; % [hr]
dt_prime = 1.5; % [hr]
t=0.5:dt:6.0;
u_0_5 = [404 1079 2343 2506 1460 453 381 274 173 0 0 0];
M = zeros(length(t),length(t));
for i=1:length(t)
u_temp = u_0_5(1:end-i+1);
M(i:12,i) = u_temp(1:end)';
end
s_hydo = dt*sum(M,2);
s_hydro_lag = [0;0;0;s_hydo(1:end-3)];
u_1_5 = 1/dt_prime*(s_hydo-s_hydro_lag);

2-3 Synthetic Unit Hydrographs

Thus far, the methods we have described only apply for the basin with gage measurements of the streamflow, however, many basins are ungauged. Therefore, engineers have analyzed many unit hydrographs from gauged basins to obtain a universal functional representation for the unit hydrograph that can be generalized to ungauges basins. This representation of the unit hydrograph is often called Synthetic Unit Hydrograph (SUH) that can be applied to ungauged basins. These synthetic methods typically scale a general UH shape or calculate key ordinates based on easily determined characteristics of watershed and excess precipitation.
There are many of these methods, such as Snyder's or Clark's SUH; however, we will cover the very popular NRCS (SCS) dimensionless synthetic hydrograph method. The dimensionless SCS hydrograph and its triangular representation are shown in Figure 9 and tabulated in Table shown as Figure 11.
Figure 9: Dimensionless curvilinear and equivalent triangular SCS unit hydrographs (USDA, 1986)

Input parameters:

2-3-1 Concentration time

As we previously discussed the time of concentration () for a watershed is the time for a raindrop to travel from the farthest point in the watershed to the outlet. The SCS recommends two methods for calculation of :
where is the average slope of the watershed in %, L denotes the hydraulic length of the main channel in feet, and [inches] represents the potential maximum retention that is obtained from the curve numbers (CN). The SCS method suggests the following relationship for the concentration time
where and V are hydraulic length and velocity in [ft] and [ft/s], respectively. The velocity of the upland method can be obtained from Figure 8 based on the average slope of the watershed.
Clearly, if we discretize the the stream into k segments, the concentration time is the sum of travel time for different stream segments as we discussed before,
.
Figure 10: Velocities for upland method for estimation of (USDA, 1986).

2-3-2 Time to Peak

Time to peak () is the time from the beginning of the rainfall to the time of the peak discharge, which can be expressed as a function of the time lag as follows:
where is in , is the duration of the rainfall excess in , and is the lag time in .
When is not known, the SCS method suggests
and because , one can obtain the following relationship for approximating the time to peak:

2-3-3 Peak Discharge

From continuity equation, we know that the area under the unit hydrograph is equal to the volume of the excess unit rainfall - example 1 [in] all over the watershed. Here, the unit of in SCS method is the unit of discharge per depth of rainfall excess such as or .
Therefore, base on conservation of mass of the depth of the rainfall and the triangular shape of the proposed SCS unit hydrograph, we have
which can be rearranged as
For the SI unit, in the above equation the area should be in and time is in second to obtain the discharge in . To convert the above formula to be used for watershed area in square kilometers and hourly time scale, we need to to multiply the above equation by .
Moreover, the SCR method suggests and thus
where . The value in English unit is 483.4, where A is in square miles.
After obtaining the and , the curvilinear coordinates of the SUH can be obtained from the table in the next slide.
Figure 11: Ratios for dimensionless SCS synthetic unit hydrograph and mass curve.

2-3-4 Construction of SCS synthetic UH

  1. Compute the time of concentration () using the lag method or the upland-velocity method.
  2. Compute the time to peak and then the peak discharge .
  3. Compute the the base time of the SUH as , which is for the triangular hydrograph and for the curvilinear one and the recession time .
  4. Compute the SUH ordinates and plot them. For the the triangular SUH, the ordinates can be directly computed while for the curvilinear one the dimensionless ratios shall be used.
Figure 12: The SCS unit hydrographs (USDA, 1986)
---------------------------------------------------------------
Example 7.2: Construct a 10-minute unit hydrograph for a basin with an area of 3.0 and concentration time [hr].
[min] = 0.166 [hr] (duration of excess rainfall)
[hr] (lag time)
[hr] (time to peak)
(peak flow discharge)
Figure 13: The calculated SCS SUH.

4- Flow Routing

Flow routing is the procedure of determining the downstream hydrograph from hydrographs at one or more points upstream.
As water enters a channel, either from hillslopes or upstream channels, the flow velocity shifts the location of the peak flow and surface roughness and storages along the water path attenuate the magnitude of the peak. The flow routing is referred to a class of methods that calculates the redistribution (widening) and translation of upstream hydrograph as a function watershed flow characteristics.
Figure 14: Left: Schematic of redistribution and translation (Chow et al. 1988). Middle: The storage accumulation and release during the flow routing. Right: Conceptualization of channels as storage reservoirs (COMET Program)
There are two main types of routing methods:
1. Lumped (Hydrologic) Routing: Flow is calculated at a particular point as a function of time alone.
2. Distributed (Hydraulic) Routing: Flow is calculated as a function of time and space.
For this course, we focus on the simpler lumped modeling.
Within lumped modeling, there are multiple methods depending on what system we are trying to model. We are going to focus on two common methods for flow routing in reservoirs (level-pool method) and stream channels (Muskingum method).
Derivation of routing methods begins with the continuity equation:
In routing, is the known inflow hydrograph, while and are unknown. Just as in our discussion of linear reservoirs, we need a storage-discharge relationship, .
There are two types of storage-discharge relationships:
Figure 15: Invariable reservoir (left) that is wide and deep compared to its length. The flow velocity is low and the water surface is almost horizontal. A variable reservoir (right) is a long, narrow and shallow reservoir with high flow velocity such that the water surface profile is markedly curved due to the backwater effect.
Figure 17: Photos of a reservoir, dam and a pond with hydraulic structures as flow controls.

4-1 Level Pool Routing

Level pool routing typically applies to man-made hydraulic structures like reservoirs, dams and stormwater ponds that are common in many watersheds. The key assumption is a horizontal water surface and negligible velocity in the reservoir. This allows an invariable relationship between storage (water surface height) and discharge, because any change in storage results in a uniform change in height of water surface.
Beginning with the continuity equation , we have , one can discretize it as follows assuming that the changes in inflow and outflow are linear over a small timestep :
Figure 17: Discretization of the reservoir continuity equation over a small time interval.
The inflow and are known at all times while initial storage and outflow are specified as initial conditions, so we can re-arrange the equation into knowns and unknowns as follows:
To use this equation, we need the elevation-storage function and elevation-discharge relationships of the reservoir to create a storage-outflow function that allows obtaining the outflow at time step from the left-hand side of the above equation iteratively.
For example, for an uncontrolled ogee crest spillway, the outflow as a function of total water head on the crest is:
is the effective length of the crest
is the discharge coefficient
is discharge in [cfs].
Figure 18: Discharge relationship for vertical-faced ogee crest (US bureau of reclamation 1987). The figure is for English system and is the water surface upstream from weir drawdown.
For a general reservoir, once the storage-outflow function is calculated, as schematically shown in Figure 19, one can calculate the RHS of Eq.(1) and compute the LHS.
Figure 19: Change of storage during a routing period and development of the storage-outflow function for level pool routing on the basis of storage-elevation-outflow curves (Chow et al. 1988).
---------------------------------------------------------------
Example 7.3: Inflow hydrograph is specified at min. Horizontal area of the reservoir is acre (43,560 ) and . Given the inflow hydrograph, compute the outflow hydrograph.
Figure 20: The left four columns represent the problem inputs and construction of the storage-outflow curve . The last six columns represent the inflow and the routed outflow. An excel worksheet is posted in the class website.
Solution:
[cfs]
Now we have to use the curve to obtain from :
[cfs] (interpolation)
[cfs]
[cfs] (interpolation)
- Computation continues for further time steps ⋯
The whole pocess can be formulated as follows:
Figure 21: The storage-outflow curve (left) and inflow and routed outflow hydrographs for the level-pool routing example.

4-2 Hydrologic River Routing

In order to calculate the routing for streams in a watershed, we can no longer assume that the water surface is level, since there is not negligible velocity in a stream and there are backwater effects. Therefore, we need to define a different storage-outflow equation, .
The Muskingum method assumes that storage in a channel is the sum of a prism and a wedge storage.
This wedge storage represents a flood wave as it propagates through a stream channel as shown in Figure 22.
Figure 22: Schematics of wedge storage and flood wave propagation assumed in Muskingum method (Credit: The COMET Program)
We can discretize this storage equation as follows for time step jth
and thus:
Recall, we previously defined the change in storage as follows:
Combining equations 2 and 3, one can obtain the following flood routing equation:, known as the Muskingum stream routing formulation:
where .
Typically, the streams is broken up into N subreaches when performing this calculation. The U.S. Army Corps of Engineers recommends the following condition to ensure solution stability when we discretize the river reach:
---------------------------------------------------------------
Example 7.4:The inflow hydrograph to a river reach is given below. Determine the outflow hydrograph of this reach, if [hr], and [hr]. The initial outflow is 85 .
Figure 23: The results of the Muskingum routing.